METHODS AND APPARATUSES OF ESTIMATING THE POSITION OF A 



MOBILE USER IN A SYSTEM OF SATELLITE DIFFERENTIAL NAVIGATION 



FIELD OF THE INVENTION 

5 The present invention relates to methods of information processing in satellite 

navigation systems with differential positioning of a mobile user. 



BACKGROUND OF THE INVENTION 

Satellite navigation systems, such as GPS (USA) and GLONASS (Russia), are 
10 intended for high accuracy self-positioning of different users possessing special 
navigation receivers. A navigation receiver receives and processes radio signals 
broadcasted by satellites located within line-of-sight distance. The satellite signals 
comprise carrier signals that are modulated by pseudo-random binary codes, which are 
then used to measure the delay relative to local reference clock or oscillator. These 
15 measurements enable one to determine the so-called pseudo-ranges between the 
receiver and the satellites. The pseudo-ranges are different from true ranges (D, 
distances) between the receiver and the satellites due to variations in the time scales of the 
satellites and receiver and various noise sources. To produce these time scales, each 
satellite has its own on-board atomic clock, and the receiver has its own on-board clock, 
20 which usually comprises a quartz crystal. If the number of satellites is large enough 

(more than four), then the measured pseudo-ranges can be processed to determine the user 
location (e.g., X, Y, and Z coordinates) and to reconcile the variations in the time scales. 
Finding the user location by this process is often referred to as solving a navigational 
problem or task. 

25 The necessity to guarantee the solution of navigational tasks with accuracy better 

than 10 meters, and the desire to raise the stability and reliability of measurements, have 
led to the development of the mode of "differential navigation ranging," also called 
"differential navigation" (DN). In the DN mode, the task of finding the user position is 
performed relative to a Base station (Base), the coordinates of which are known with the 

30 high accuracy and precision. The Base station has a navigation receiver that receives the 
signals of the satellites and processes them to generate measurements. The results of 
these measurements enable one to calculate corrections, which are then transmitted to the 
user that also uses a navigation receiver. By using these corrections, the user obtains the 



ability to compensate for the major part of the strongly correlated errors in the measured 
pseudo-ranges, and to substantially improve the accuracy of his or her positioning. 

Usually, the Base station is immobile during measurements. The user may be 
either immobile or mobile. We will call such a user "the Rover." The location 
5 coordinates of a moving Rover are continuously changing, and should be referenced to a 
time scale. 

Depending on the navigational tasks to be solved, different modes of operation 

may be used in the DN mode. They differ in the way in which the measurement results 

are transmitted from the Base to the Rover. In the Post-processing (PP) mode, these 
10 results are transmitted as digital recordings and go to the user after all the measurements 

have been finished. In the PP mode, the user reconstructs his or her location for definite 

time moments in the past. 

Another mode is the Real-Time Processing (RTP) mode, and it provides for the 

positioning of the Rover receiver just during the measurements. The RTP mode uses a 
1 5 communication link (usually it is a radio communication link), through which all the 

necessary information is transmitted from the Base to the Rover receiver in digital form. 
Further improvement of accuracy of differential navigation may be reached by 

supplementing the measurements of the pseudoranges with the measurements of the 

phases of the satellite carrier signals. If one measures the carrier phase of the signal 
20 received from a satellite in the Base receiver and compares it with the carrier phase of the 

same satellite measured in the Rover receiver, one can obtain measurement accuracy to 

within several percent of the carrier's wavelength, i.e., to within several centimeters. 

The practical implementation of those advantages, which might be guaranteed by 

the measurement of the carrier phases, runs into the problem of there being ambiguities in 
25 the phase measurements. 

The ambiguities are caused by two factors. First, the difference of distances AD 

from any satellite to the Base and Rover is much greater than the carrier's wavelength X. 

Therefore, the difference in the phase delays of a carrier signal Aqp^ADIX received by the 

Base and Rover receivers exceeds several cycles. Second, it is not possible to measure 
30 the integer number of cycles in Aq) from the incoming satellite signals; one can only 

measure the fractional part of A(p. Therefore, it is necessary to determine the integer part 

of A(p, which is called the "ambiguity". 
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More precisely, we need to determine the set of all such integer parts for all the 
satellites being tracked, one integer part for each satellite. One has to determine this set 
along with other imknown values, which include the Rover's coordinates and the 
variations in the time scales. 
5 In a general way, the task of generating highly-accurate navigation measurements 

is formulated as follows: one determines the state vector of a system, with the vector 

containing unknown components. Those include three Rover coordinates (usually 

along Cartesian axes X, Y, Z) in a given coordinate system (sometimes time derivatives 
of coordinates are added too); the variations of the time scales which is caused by the 

10 phase drift of the local main reference oscillator; and n integer imknown values associated 
with the ambiguities of the phase measurements of the carrier frequencies. The value of n 
is determined by the nvunber of different carrier signals being processed, and accordingly 
coincides with the number of satellite channels actively functioning in the receiver. At 
least one satellite channel is used for each satellite whose broadcast signals are being 

1 5 received and processed by the receiver. Some satellites broadcast more than one code- 
modulated carrier signal, such as a GPS satellite that broadcasts a carrier in the Li 
frequency band and a carrier in the L2 frequency band. If the receiver processes the 
carrier signals in both of the Li and L2 bands, the number of satellite channels (n) 
increases correspondingly. 

20 Two sets of navigation parameters are measured by the Base and Rover receivers, 

respectively, and are used to determine the unknown state vector. Each set of parameters 
includes the pseudo-range of each satellite to the receiver, and the full (complete) phase of 
each satellite carrier signal, the latter of which may contain ambiguities. Each pseudo- 
range is obtained by measuring the time delay of a code modulation signal of the 

25 corresponding satellite. The code modulation signal is tracked by a delay-lock loop 
(DLL) circuit in each satellite-tracking channel. The ftiU phase of a satellite's carrier 
signal is tracked by phase counter (as described below) with input from a phase-lock-loop 
(PLL) in the corresponding satellite tracking chaimel (an example of which is described 
below in greater detail). An observation vector is generated as the collection of the 

30 measured navigation parameters for specific (definite) moments of time. 

The relationship between the state vector and the observation vector is defined by 
a well-known system of navigation equations. Given an observation vector, the system of 
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equations may be solved to find the state vector if the number of equations equals or 
exceeds the number of unknowns in the state vector. In the latter case, conventional 
statistical methods are used to solve the system: the least-squares method, the method of 
dynamic Kalman filtering, and various modifications of these methods. 
5 Practical implementations of these methods in digital form may vary widely. In 

implementing or developing such a method on a processor, one usually must find a 
compromise between the accuracy of the results and speed of obtaining results for a given 
amount of processor capability, while not exceeding a certain amount of loading on the 
processor. The present invention is directed to novel methods and apparatuses for 
10 accelerating the obtaining of reliable estimates for the integer ambiguities at an acceptable 
processor load. 

More particularly, the present invention is directed to novel methods and 
apparatuses for more quickly obtaining such estimates in floating-point form (non-integer 
form) which are close to the integer values. With these floating-point forms, which we 
1 5 call floating ambiguities, conventional methods may be used to derive the corresponding 
integer ambiguities. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is schmatic diagram of an exemplary receiver which may be used in 
practicing the present invention. 
20 FIG. 2 is a perspective view of a rover station and a base station in a first 

exemplary coinfiguration according to the present invention. 

FIG. 3 is a perspective view of a rover station and a base station in a second 
exemplary coinfiguration according to the present invention. 

FIG. 4 is a flow diagram illustrating a general set of exemplary embodiments 
25 according to the present invention. 

FIG. 5 is a schematic diagram of an exemplary apparatus according to the present 
invention. 

FIG. 6 is a diagram of an exemplary computer program product according to the 
present invention. 

30 FIG. 7 is a flow diagram illustrating another general set of exemplary 

embodiments according to the present invention. 

FIG. 8 is a diagram of another exemplary computer program product according to 
the present invention. 
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SUMMARY OF THE INVENTION 

Broadly stated, the present invention encompasses methods and apparatuses for 
estimating the floating ambiguities associated with the measurement of the carrier signals 
5 of a plurality of global positioning satellites, such that the floating ambiguities are 
preferably consist for a plurality of different time moments. 

The floating ambiguities are associated with a set of phase measurements of a 
plurality n of satellite carrier signals made by a first navigation receiver (B) and a second 
navigation receiver (R) separated by a distance, wherein a baseline vector {x^,y^,z^) relates 
10 the position of the second receiver to the first receiver. Each satellite carrier signal is 
transmitted by a satellite and has a wavelength, and each receiver has a time clock for 
referencing its measurements. Any difference between the time clocks may be 
represented by an offset. Methods and apparatuses according to the present invention 

receive, for a plurality of two or more time moments 7, the following inputs: a vector yf 

1 5 representative of a plurality of pseudo-ranges measured by the first navigation receiver 
(B) and corresponding to the plurality of satellite carrier signals, a vector ^ 
representative of a plurality of pseudo-ranges measured by the second navigation receiver 
(R) and corresponding to the plurality of satellite carrier signals, a vector Dj 
representative of a plurality of estimated distances between the satellites and the first 

20 navigation receiver (B), a vector Dj representative of a plurality of estimated distances 
between the satellites and the second navigation receiver (R), a vector q)j representative 
of a plurality of full phase measurements of the satellite carrier signals measured by the 
first navigation receiver (B), a vector qjj representative of a plurality of full phase 
measurements of the satellite carrier signals measxired by the second navigation receiver 

25 (R), and a geometric Jacobian matrix whose matrix elements are representative of the 
changes in the distances between the satellites and one of the receivers that would be 
caused by changes in that receiver's position and time clock offset. (As used herein, the 
term "representative of," as used for example used when indicating that a first entity is 
representative of a second entity, includes cases where the first entity is equal to the 

30 second entity, where the first entity is proportional to the second entity, and where the 
first entity is otherwise related to the second entity.) 
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The present invention may be practiced in real time, where estimates of the 
floating ambiguities are generated as the above satellite information is being received. 
The present invention may also be practiced in post-processing mode, where the floating 
ambiguities are estimated after all the above satellite information has been received. In 
5 the latter case, block processing according to the present invention may be done. 

Preferred methods and apparatuses according to the present invention generate, for 

each time moment y, a vector A}^- of a plurality of range residuals of pseudo-range 
measurements made by the first and second navigation receivers in the form of: 

^rj= {'^-Tf)-{Dj^-Dfy, 

10 and also generate, for each time moment y, a vector ^ of a plurality of phase residuals 
of full phase measurements made by the first and second navigation receivers in the form 

of: Aq>i-iq>^^ (p/) - A"' •(/>/ - />/ ), 

, -1 

where A is a diagonal matrix comprising the inverse wavelengths of the satellites, hi 
the case of real-time processing, an LU-factorization of a matrix Mi, or a matrix inverse 
15 of matrix Mi, is generated for a first time moment (denoted as j = 1), with the matrix Mi 

being a function of at least A and Hj^. Also for this initial time moment, an initial 
vector Ni of floating ambiguities is generated as a function of at least A//, A <pj, and the 
LU-factorization of matrix Mi or the matrix inverse of matrix Mi. For an additional time 
moment j, an LU-factorization of a matrix My, or a matrix inverse of matrix My, is 

20 generated, with the matrix My being a function of at least A and Also for an 

additional time moment j, a vector Ny of estimated floating ambiguities is generated as a 
function of at least A^-, A g)j , and the LU-factorization or matrix My or the matrix inverse 
of matrix My. Exemplary forms of matrices My and vectors Ny are provided below. In this 
manner, a set of successively more accurate estimates of the floating ambiguities are 

25 generated in real-time with the vectors Ny. This method, of course, may also be practiced 
in a post-processing environment, where the data has been previously recorded and then 
processed according to the above steps. 

In the post-processing environment, the following block processing approach 
according to the present invention may be practiced. As with above-described 

30 embodiments for real-time processing, the vectors of pseudo-range residuals Ay^ and 



6 



vectors of phase residuals 2I ^Jj^ , k=l,. . .7, are generated. Thereafter, a general matrix M 

is generated from the data, with M being a function of at least A and H^^ for index k 
of J3jt^ covering at least two of the time moments 7, and an LU- factorization of matrix M 
or a matrix inverse of matrix M is generated. Thereafter, a vector N of estimated floating 
S ambiguities is generated as a ftmction of at least the set of range residuals A^, the set of 

phase residuals A q)^^ and the LU- factorization of matrix M or the matrix inverse of 
matrix M. 

As an advantage of the present invention, the nature of matrix M, as described in 
greater detail below, enables a compact way of accumulating the measured data in order 
10 to resolve the floating ambiguities. As a fiirther advantage, forms of matrix M provide a 
stable manner of factorizing the matrix using previous information. In preferred 
embodiments, matrix M is substantially positive definite, and also preferably symmetric. 

Accordingly, it is an objective of the present invention to improve the stability of 
generating estimates of the floating ambiguities, and a fiirther objective to reduce the 
15 amount of computations required to generate estimates of the floating ambiguities. 

This and other advantages and objectives of the present invention will become 
apparent to those of ordinary skill in the art in view of the following description. 

DETAILED DESCRIPTION OF THE PRESENT INVENTION 

20 

NOMENCLATURE 

All vectors presented herein are in column form. The transpose of a vector or matrix is 
indicated by the conventional superscript affixed to the end of the name of the vector 
or matrix. The inverse of a matrix is indicated by the conventional superscript 
25 affixed to the end of the name of the matrix. For the convenience of the reader, we 
provide a summary notation symbols below: 

y " denotes a pseudorange measurement, also called a code measurement; a pseudorange 
is the approximate distance between a satellite and a receiver is equal to the speed of light 
C multiplied by the time it takes the satellite signal to reach the receiver. 
30 (p - denotes a phase measurement. 

n ' denotes the number of satellites fi"om which measurements are collected 
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- is a superscript index that is attached to a measurement quantity (e.g,, y , <p)or 

other data to identify the satellite to which the quantity corresponds. 

j and k - are subscript indices that are attached to a measurement quantity (e.g., y y <p) 

or other data to identify the moment in time {e.g., epoch) to which the measure the 
5 quantity corresponds, 
c - speed of light in air. 

X - denotes a wavelength of a carrier signal emitted by the S-th satellite. The 

wavelength may denote either the LI -band carrier or the L2-band carrier, as indicated by 
the context of the discussion. 

10 A - denotes a diagonal matrix of wavelengths: A = diag{?} , • • • , 2" ) . 

O/^ orO/x;;^ - denotes an Ixm zero matrix (all matrix elements are zero). The indices 

/ and m take on the appropriate values in the context that the matrix is used. In other 
words, specific numerical values or other indices may appear in the places of I and m. 
1/ or I/x/- denotes a square identity matrix of size Ixl . The index / takes on the 
15 appropriate value in the context that the square identity matrix is used, hi other words, 
specific numerical values or other indices may appear in place of/. 
|i - is a vector of observable equations, preferably linearized observable equations. 

fi = (a;^^ , • • ■ , A/ " , A<p^ , • • • , Aip^ y , where each element of ji is associated with a 
corresponding satellite, where the superscript indices affixed to the elements identify the 
20 sateUites. 

- denotes the observation matrix H associates with the pseudo-ranges y. In the art, 

is also called the Jacobian matrix, the Jacobi matrix, the geometric function matrix, the 
Jacobian geometric matrix, the matrix of directional cosines, and the directional cosine 

matrix. H{ has dimensions of nx4. 
having dimensions of 2n x 4 . 
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25 Qk - denotes a compoimd matrix formed by and A as follows: Q. = 



G - denotes a compound matrix G = 



O 



fomied by the zero matrix 0„ „ and the 



identity matrix I„ , having dimensions of 2« x w . The matrix G has many purposes. One 
purpose is to select the lower-right nxn sub-matrix of a larger 2n x 2n matrix as follows: 

Cii I Ci2 



G'^ C G = C22, where C = 



.C21 I C22 . 



5 ' denotes a compound matrix = [Qj I ^] , whichis a2wx(4 + w) matrix. 



II • 1^-1 denotes F"^ - weighted norm, ||x||^.i = X^F"^X. 



Brief Background on the Structure of the Satellite Signals 

10 Before describing the present invention, we briefly describe the structure of the 

satellite signals and of a typical receiver suitable for differential navigation applications. 
Each of the satellites radiates signals in two frequency bands: the LI band and the L2 
band. Two carrier signals are simultaneously transmitted in the LI -band; both carrier 
signals have the same frequency, but are shifted in phase by ^/2 (90*^). The first LI carrier 

15 signal is modulated by the clear acquisition C/A-code signal and the second LI carrier 

signal is modulated by the precision P-code signal. One carrier signal is transmitted in the 
L2 band, and uses a different frequency than the LI carrier signals. The L2 carrier signal 
• is modulated by the same P-code signal used to modulate the second LI carrier signal. 
These carrier frequencies are between 1 GHz and 2 GHz in value. Each C/A-code signal 

20 and P-code signal comprises a repeating sequence of segments, or "chips", where each 

chip is of a predetermined time period (A) and has a pre-selected value, which is either +1 
or -1 . The segment values follow a pseudo-random pattem, and thus the C/A-codes and 
the P-codes are called pseudo-random code signals, or PR-code signals. Additionally, 
before each C/A-code signal and P-code signal is modulated onto its respective carrier 

25 signal, each code signal is modulated by a low frequency (50 Hz) information signal (so- 
called information symbols). 

The approximate distance between a satellite and a receiver is equal to the speed of light c 
multiplied by the transmission time it takes the satellite signal to reach the receiver. This 
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approximate distance is called the pseudorange 7, and it can be corrected for certain errors 
to find a corrected distance D between the satellite and the receiver. There is a 
pseudorage between each visible satellite and the receiver. The transmission time fi-om 
satellite to receiver is measured with the aid of clocks in the receiver and the satellite, and 
5 with the aid of several time scales (/.e., timing marks) present within the received satellite 
signal. The clocks in the satellites and the receiver are set to substantially the same time, 
but it is assumed that the receiver clock has a time offset T because the receiver clock is 
based upon a quartz-crystal whereas each satellite clock is based upon a more accurate 
atomic reference clock. The receiver has the orbital pattems of the satellites stored in a 

10 memory, and it can determine the orbital position of the satellites based on the time of its 
clock. The receiver reads the timing marks on the satellite's signal, and compares them 
against it own clock to determine the transmission time from satellite to receiver. The 
satellite's low-frequency (50Hz) information signal provides the least precise timing 
information, the C/A-code signal provides the next most precise timing information, and 

15 the P-code signal provides the most precise timing information. The pseudorange is 
determined from the low-frequency information signal and the C/A-code signal for 
civilian users and some military users, and is determined from the low-frequency 
information signal and the P-code signal for most military users. Accurate use of the P- 
code signal requires knowledge of a certain code signal that is only known to military 

20 users. Precision better than that provided by the P-code signal can be obtained by 

measuring the phase of the satellite carrier signal in a differential navigation mode using 
two receivers. 

Referring to FIG. 1, a typical receiver for differential navigation applications has 
an antenna, an amplifier/filter imit that receives the antenna's output, a frequency down- 

25 conversion unit that receives the output of the amplifier/filter imit, and several individual 
tracking channels of the coherent type, each of which receives the down-converted 
satellite signals. The receiver also comprises a receiver clock, which provides base 
timing signals to the local oscillator of the down-conversion unit, and to components 
within each individual tracking channel. The down-conversion unit provides down- 

30 converted and quantized versions of the satellite signals to the channels, with the down- 
converted signals having frequencies generally in the range of 10 MHz to 20 MHz. Each 
channel tracks a selected one of the satellite signals. Each tracking channel measures the 
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delay of one PR-code signal within the sateUite signal (e,g,, C/A-code or P-code signal), 
and also the phase of the down-converter version of the satellite's carrier signal. A typical 
tracking channel comprises a Delay-Lock Loop (DLL) circuit for tracking the delay of the 
PR-code, a Phase-Lock Loop (PLL) circuit for tracking the phase of the satellite's carrier 
5 signal, and three correlators which generate the input signals for the DLL and PLL 
circuits. 

Referring to FIG. 1, the DLL circuit has a reference code generator that generates 
a set of reference code signals, each of which tracks the PR-code of the satellite signal, 
and each of which is provided as an input to a respective correlator. Each correlator 

10 output is representative of the degree to which the reference code signals are tracking the 
satellite code signal (i.e., the amount by which the reference signals are advanced or 
retarded in time with respect to the satellite signal). The DLL circuit also has a DLL 
discriminator and a DLL filter. The DLL discriminator receives inputs from two 
correlators, a DLL correlator that generates a base signal for controlling the DLL circuit 

15 and a main correlator that generates a signal useful for normalizing the base signal from 
the DLL correlator. The DLL discriminator generates an error control signal from these 
inputs; this error signal is filtered by the DLL filter before being provided to the DLL 
reference code generator as a control signal. The value of the DLL error signal varies by 
the amount that the reference code signals of the DLL generator are delayed or advanced 

20 with respect to the satellite code signal being tracked, and causes the code generator to 
adjust the timing of the reference code signals to better match the satellite's code signal. 
In this manner, tracking is achieved. The pseudorange 7may be generated by methods 
known to the art from the receiver's clock, the 50 Hz information signal, and any one of 
the reference code signals generated by the DLL. This is indicated in the "')5(Tj)" box. 

25 In a similar manner, the PLL has a reference carrier generator that generates a 

reference carrier signal that tracks the down-converter version of the satellite's carrier 
signal. We denote the frequency of the reference carrier signal as /nco since the reference 
carrier frequency is often generated by a numerically-controlled oscillator (NCO) within 
the reference carrier generator. Referring to FIG. 1, the PLL circuit also has a PLL 

30 discriminator and a PLL filter. The PLL discriminator receives inputs from two 

correlators, a PLL correlator that generates a base signal for controlling the PLL circuit, 
and the main correlator that generates a signal useful for normalizing the base signal from 
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the PLL correlator. Each correlator output is representative of the degree to which the 
reference carrier signal is tracking the satellite's carrier signal (i.e., the amount by which 
the reference carirer signal is advanced or retarded in time with respect to the satellite's 
carrier signal). The PLL discriminator generates an error control signal from these inputs; 
this error control signal is filtered by the PLL filter before being provided to the PLL 
reference carrier generator as a control signal. The value of the PLL error signal varies by 
the amount that the reference carrier signal of the PLL reference carrier generator is 
delayed or advanced with respect to the satellite carrier signal being tracked, and causes 
the reference carriergenerator to adjust the timing of the reference carrier signal to better 
match the satellite's carrier signal. The reference carrier generator provides a phase signal 



S^jNCoiTj) representative of the phase (integrated frequency) of the PLL's reference 

oscillator. This signal may be combined with other information, as described below, to 
provide a signal representative of the full phase of the satellite signal (but possibly with 

ambiguities in it). This is indicated in the "</li(Tj)" box. 

Finally, each individual tracking channel usually comprises a search system which 
initially varies the reference signals to bring the PLL and DLL circuits into lock mode. 
The construction of this and the other above components is well known to the art, and a 
fiirther detailed description thereof is not needed in order for one of ordinary skill in the 
art to make and use the present invention. 

Brief Background on the Navigation Parameters. 

Because of the time offset rand other error effects, the pseudorange ybetween a 
satellite and a receiver is usually not the actual distance between the satellite and the 
receiver. By looking at the pseudoranges of four or more satellites, there are well-known 
methods of estimating the time offset T of the receiver and of accoimting for some of the 
error effects to generate the computed distance D between the satellites and the receiver. 
The receiver's position may then be computed. However, because of various sources of 
noise and the relatively low resolution of the pseudo-random code signal, the true 
distances (i.e., true ranges), and receiver's position coordinates will not be exactly known, 
and will have errors. 

In theory, more precise values for the receiver's position and clock offset could be 
obtained by measuring the number of carrier cycles that occur between each satellite and 




12 



the receiver. The phase of the carrier of the sateUite signal as transmitted by a satellite 
can be expressed as: 

<l>\t) = ^o -^]fs'dt 
0 

= ^0 + fs't (2) 

where is an initial phase value, fs is the satellite carrier frequency, and t is the 
satellite's time. Because the satellite carrier frequency is generated from a very precise 
time base, we may assume that is a constant and not time dependent, and we may 
replace the above integral with fs^t, as we have shown in the second line in the above 
equation. The phase of this signal when it reaches the receiver's antenna over the range 
distance D{t) is denoted as <f>^AOX its value would be: 
</>^^ (0 = <l>\t - D(t)/c - T^um 

^^O'^fsit- D{t)lc - r^rJi) ) 

-^^0 +fst -fs'D(t)/c -fsT,Ut\ (3) 
where c is the speed of light, and where TatJJ) is a delay due to anomalous atmospheric 
effects which occur mostly in the upper atmosphere. The number of cycles fs'TATiAf) due 
to the atmospheric effects cannot be predicted or determined to within an accuracy of one 
carrier cycle by a stand-alone receiver (/.e., cannot be determined by absolute 
positioning). However, the atmospheric effect can be substantially eliminated in a 
differential GPS mode where the phase of the satellite is measured at a rover station and a 
base station, 4^ A,Rif) and <P^aMO respectively, and then subtracted from one another. 
Over the short baseline between the rover and base stations, the atmospheric delay T^rjiXO 
in both of these phases is equal for practical purposes, and the difference in phases is: 
<I>\r(0 - <I>^aAO =fs'DR(t)/c - fsDB{t)lc. (4) 
The terms $ q and fs't have also been cancelled from the difference. In FIG. 2, 
we show a base station, a rover station, and four satellites, in schematic form. The wave 
fronts of the satellite carrier signal are nearly planar by the time they reach the base and 
rover stations since the satellites are at least 22,000 km away. These wave fronts are 
illustrated in FIG. 2 for the first and fourth sateUites (SI and S4). Assuming that a wave 
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front of the satellite carrier reaches one of the antennas first (either rover or base), the 
above difference is the number of additional wave fronts that the satellite carrier must 
travel before it reaches the second antenna (either base or rover). We have illustrated this 
for the first satellite SI as fsi'DR{i)/c - fsi'DB(t)/c, Using the angle between the base- 
5 rover baseline vector and the Une-of-sight vector to the satellite from either of the 
stations, the number of carrier cycles that lie along the baseline can be determined by 
trigonometry, and thus the baseline can, in theory, be accurately determined from the 
differential phase. 

However, the task of measuring carrier phase is not as easy as it appears. In 

10 practice, we must use non-ideal receivers to measure the phases 4^A,Rif) and 4Pa,b{^\ with 
each receiver having a different clock offset with respect to the GPS time, and with each 
receiver having phase errors occurring during the measurement process, hi addition, at 
the present time, it is not practical to individually count the carrier cycles as they are 
received by the receiver's antenna since the frequency of the carrier signal is over 1 GHz. 

1 5 However, the PLL loop can easily track the Doppler-shift frequency fo of the carrier 
signal, which is in the kHz range. With a few assumptions, the phases (t>^A,R(0 and 
4^A,B{t) can be related to their respective Doppler-shift frequencies. As is known in the 
art, the satellite transmits at a fixed frequency of fs^ but the relative motion between the 
satellite and receiver causes the frequency seen by the receiver to be slightly shifted from 

20 the value of^ by the Doppler frequency fr,. We may write the frequency seen by the 

receiver's antenna ^sfs-^fo^ where fo has a positive value when the distance between the 
satellite and receiver's antenna is shrinking, and a negative value when the distance is 
increasing. Each receiver can then assume that the received phase is proportional to the 
predictable amount of fs't, minus the amount fo't due to the Doppler-shift. The Doppler 

25 amount fo't is subtracted from fs't because the Doppler frequency increases as the 

distance between the satellite and receiver's antenna decreases. The predictable amount 
fs't will be the same for each receiver, but the Doppler frequencies will usually be 
different. 

As previously mentioned with reference to FIG. 1, the reference oscillator (e.g., 
30 NCO) of the PLL circuit tracks the frequency of a selected one of the down-converted 
satellite signals. As a result, it inherently tracks the Doppler-shift frequency of the 
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satellite's carrier signal. Before being provided to the PLL circuit, the carrier signal is 
down-converted from the GHz range by a local oscillator having a frequency/t,. The 
frequency seen by the PLL circuit is (/^ + f^fi^ which can be rearranged as: {fs -fi) + 
fo- The quantity (fs -fi) is known as the pedestal frequency fp, which is typically around 
5 10 MHz to 20 MHz. The PLL's reference oscillator tracks the down-converted frequency 
oifp+fo^ We would like to integrate the frequency of the reference oscillator to obtain a 
phase observable which is proportional to - fo^t. Starting at a time moment Tp when the 
PLL circuit initially locks onto the down-converted satellite signal, with the time moment 
Tp being measured by the receiver's clock, we will generate a phase observable ^jiTj) at 

10 discrete time moments Tj of the receiver's clock T'as follows: 

<t>,iTj) -fp.nom\Tj - Tp) - <j>jNC0(Tj) (5) 

where fp^nom is the nominal value of the pedestal frequency, and where <l>jNCoiTj) is the 
phase (integrated frequency) of the PLL's reference oscillator {e,g.^ NCO). The time 
moments Tj are spaced apart from each other by a time interval A7y, as measured by the 

15 receiver's clock, and may be express asTj — j'^Tj , where j is an integer. (f>j(Tj) is in 
units of cycles, and is proportional to the negative of the integrated Doppler-shift 
frequency. This is because <l>jNCo{Tj) changes value in proportional to the quantity (fp + 

fonTj-Tp). 

While it may be possible to set <l>jNCo(Tp) to any value at the initial lock moment 
20 Tp, it is preferably to set <l>jNCoiTp) to a value substantially equal to (fp,nom'Tp- 

fs' T)p /c), where Tp is measured from the start of GPS time, and where Dp is the 

approximate distance between the satellite and the receiver, as found by the pseudorange 
measured by the DLL, or as found by a single point positioning solution. This setting of 

4^jNCo{Tp) is conventional and provides values of <l>j{Tj) for the base and rover stations 

25 which are referenced from the same starting time point. 

In U.S. Patent No. 6,268,824, which is commonly assigned with the present 

application, it is shown how the observable <f>j(Tj) is related to the distance D between the 
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receiver and the satellite. We summarize the relationship here, neglecting the 
atmospheric delay TAn/jtj) since this delay will be canceled out when we take the 
difference between receivers, and refer the reader to the patent for the derivation: 

0/7}) =fsDjlc +fsiTj- Tp) + (0^,- ^1 )-N„co - (6) 

where 

• fs is the satellite carrier frequency, 

• Dj is the distance between satellite and receiver at time moment Tj, 

• c is speed of light in the atmosphere, 

• Tj is the offset between the sateUite time clock t and the receiver time clock T 
(Tj = tj -^Tj ), Tj varies with time, 

• Tp is the time moment when the PLL circuit initially locks onto the down- 
converted satellite signal, 

nom'Tp^ where $ lq is the phase offset of the receiver's local 

oscillator and fL,nom is the nominal frequency of the receiver's local oscillator, 

• ^0 IS an initial phase value of the satellite, 

• Nnco is an unknown integer number of carrier cycles due to the PLL not 
knowing which carrier cycle it first locks onto, Nnco may be positive or 
negative, and 

• ^{fij is a small tracking error due to the operation of the PLL circuit. 

As it turns out, if we chose the time increments A7y for Tj (Tj = J'ATj ) to be equal to 

ATj = 2 ms or an integer multiple of 2 ms, then the term -fs'Tp is an integer number 
which can be consolidated with the integer ambiguity -Nnco into a single ambiguity +N 

^iTj) -fs'Djlc ^fsTj + (</)J^- *o ) + A^- f0y (7) 
Thus, our Doppler phase observable <i>jiTj) has been related to the distance Dj between 
the satellite and receiver, the time offset Tj of the receiver, an integer ambiguity A^, and 
some initial phase offsets (0o- #o ) which can be readily determined. 
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We will now write equation (7) for the base and rover stations, adding superscripts 
"jB" and "i?" for the base and rover stations, and subscript "/w" to indicate the iw-th 
satellite signal and the w-th individual tracking channel. 

<t>lmaj) =fsDljc -^fs -r/ + (0^^ - $^^) + at/ - (8A) 

<l>lm{Tj)=fsDljc+fs'r^ +((l>'o''-^o)+Nj-^^'^rn (8B) 
For the differential navigation mode, the difference of these phases is formed: 
<f>j,miTj) - (j>lmiTj) = fsD^j,Jc -fsD%lc ^fs -(r/ - r/ ) 

+ (iV/ - nJ") + i(/>'o^- (i>'o'') - - ^<^J.m) (9) 

Using = (Nn^ - nJ^) to represent the difference between the ambiguities, and using 
the well-known relationship fsic =l/)wi» where >Sn is the wavelength of the satellite carrier 
signal, we have: 

^Im^Tj) - ^lm{Tj) = il/JLy(D% - D%) +fs -(if -Tf)+N'" 

+ (<f>o^ - <t>'o'') - i^Am - ^<^j,m) (10) 

The values for </)J^^ and can be readily determined. The values of and 

caimot be determined, but they have zero mean values (and should average to zero over 
time). 

Now, we relate this back to the initial objective and problem that we expressed 
with respect to FIG. 2. As indicated, we wanted to find the phase difference quantity 
(fs'DR{i)/c -fs'DB(t)/c) associated with each satellite in order to determine the number of 
carrier cycles that lie along the baseline, and thus provide a better estimate of the 
coordinates of the baseline. By examination of the phase difference quantity and 

equations (9) and (10), we seek that (fsDR{t)/c - fsDB{t)lc) is equal to ^j^rni^j) ~ 0/ 
^m(^j) y which we can measure, minus the cycle ambiguity N^y minus the error in the 
clocks/5 -(rZ-r/), minus the phase offsets (</)Jj^ - 4^o^\ ^nd minus the errors (Sif>j,m ~ 

^(f)j\m)^ the last of which can be removed by averaging. Thus, estimating the cycle 

ambiguities N will be important in using the phase information to improve the estimated 
coordinates of the baseline. 
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Our Inventions related to Ambiguity Resolution - Estimating Floating Ambiguities 

Our inventions can be applied to a number of application areas. In general 
application areas, the rover and base stations can move with respect to one another. In 
one set of specific application areas, the base and rover stations are fixed to a body or 
5 vehicle (e.g., planes and ships) and separated by a fixed distance L, as shown in FIG. 3. 
We refer to these applications as the constant distance constraint applications. The 
movement of the vehicle may be arbitrary while the distance between two antennas is held 
constant due to the rigidity of the vehicle body. This distance Lrb is known with a certain 
degree of confidence is considered as an additional constraint or measurement available at 
10 every time instant, because it depends on both Base and Rover positions. The constraint 
can be mathematically expressed as: 

(i^R -^B? ^{yR -ysf +(^7? -^bY)^^ =Lrb (11) 

where Xr, yR, and Zr, represent the rover station coordinates and where Xb, yB> Zb 
represent the base station coordinates. The information firom the constraint can be used to 
15 better estimate the Rover position relative to the Base. For this, it is preferably to form a 
cost penalty fimction mathematically equivalent to the form: 
F{xr -XB.yR -yB>ZR -Zb) = 

where, in preferred embodiments, the cost fimction F(*) will be minimized such that and 
20 its the first derivatives with respect to the rover coordinates (or alternatively the base 
coordinates) are reduced to values near zero: 

F(.).0,^»0,|^«0.^«0 

dxR dyR dzR 

Jn the preferred methods described below, the amount of weighting given to minimizing 
the cost fimction F(*) will set by a weighting parameter q (e.g., g'F(*)) , with no 
25 weighting being given when ^=0, and increasing degrees of weighting being given by 
selecting q>0. The case of setting q=0 is also equivalent to removing the constraint that 
the rover and base stations are separated by a fixed distance L. 

The ambiguity resolution task generally comprises the following three main parts: 
1 . Generation of the floating ambiguity estimations (estimating the floating 
30 ambiguities). 
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2. Generation of the integer ambiguities, and 

3. Formation of a signal of integer ambiguities resolution. 

The present invention pertains to the first part, that of generating the floating ambiguity 
estimates. The second and third parts may be accomplished by apparatuses and methods 
known in the art. The apparatuses and methods of our present invention use as an input 
data: 

(a) the pseudo-ranges and full phase measurements obtained in the Base and Rover 

receivers; 

(b) the satellites coordinates (which are highly predictable can be readily 
determined fi-om the GPS time and information provided in the 50 Hz low- 
fi-equency information signal by methods well known to the art); 

(c) the estimated coordinates of the Base and Rover stations (as measured at their 
antennas); which may be generated firom then pseudo-range data by single point 
solutions, and oftentimes the coordinates of the base station are precisely 
known; 

(d) optionally, a set of weight coefficients characterizing the accuracy of the 
measurements; and 

(e) optionally, an estimate 'of the baseline coordinates r, and the value of Lrb if 
the cost function F(*) is used. 

According to the input data, a vector of observations and a covariance matrix of 
measurements are formed. After that a state vector is generated, components of which are 
floating ambiguities, the number of which is equal to the niraiber of satellites. On the 
basis of the floating ambiguity values, a search of the integer ambiguities is performed 
with use of a least-squares method that is modified for integer estimations. 

hnprovement of the floating ambiguity estimations can take place step-by-step, 
and the probability of correct ambiguity resolution increases step-by-step as information 
is accumulated. Preferred finishing of this process is registered by appearance of the 
signal of integer ambiguities resolution, which indicates that ambiguity resolution was 
performed sufficiently safely. After that, the integer ambiguities together with other input 
data are used for accurate determination of the base line vector. The tasks of determining 
the integer ambiguities and of generating the signal of integer ambiguities resolution are 
known to the art and do not form part of our invention. These tasks, therefore, are not 
described in greater detail. 
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Our invention is suitable for both the RTP mode and for the PP mode under 
movable Rover, where the Rover coordinates may be random and independent in adjacent 
clock moments. Our invention is also suitable for both the RTP mode and for the PP 
mode when both the base and rover stations are moving, where the coordinates of the 
stations may be random, but constrained by fixed distance Lrb, in adjacent clock 
moments. 

We start by rewriting equation (10), which has the integer ambiguities in 
vector form: 

(Ihin - <l>^j^Td = A"' (jfj - Z)J) +fs (t/ - t/ ) + 

+ (<^'/-0'/)-a4-f^iX (13) 

. -1 

where A is a diagonal matrix of inverse wavelengths. We note that the components of 

(0 - <^ can be selected to be integer constants, and can therefore be incorporated 

into the integer ambiguities N, which may be called the modified integer ambiguities A^. 
Thus, we may simplify the above equation: 

<A>y) - <A/7;) = A'' •(/>* - D^) +fs -(t/ - t/ ) + iv- (r4- Uj) (14) 

After this modified N are found, the "true" N may be found by subtracting (<I>q^ q^) 
fi-om the modified A^. 

As an alternative equivalent to the above, the vector (0 ~ 0 ^?^) 
computed and carried over to the left-hand side of equation (14) and subtracted fi-om 
(<j>jiTi) - <j>jiTi)), As such, the that is determined will be the "true" N. Once the 
"true" N is found, then the underlying unknown integer number of carrier cycles Nnco 
may be found by negating the "true" N (i.e. , multipling by -1) and then subtracting/^' Tp 
(Nnco = - "true" N-fs'Tp). However, for the ultimate goal of computing more accurate 
coordinates of the baseline, the modified ambiguities in equation (14) can be used, and 
they greatly simplify the computations. Nonetheless, it will be appreciated that the 
present invention may be practiced by using the "true" N or N„co in equation (14) by 
suitable modification of the left-hand side of equation (14), and that the appended claims 
cover such practices. 
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Our approach comprises solving the above equations (14) at a plurality of time 
moments jointly to find an estimation for vector N. We cannot readily find measured 

B J7 ^ ^ 

values for vectors Tj , 7j , f , and f , and we have some errors in the range vectors 

BR B 

D j and D j. Our approach is to represent (model) the errors in the range vectors Dj and 

D j and the terms fs ' (jj -Tj ) - (f - f ) by the following error vector: A • 

Hj^\Ax, Ay, Az, c'At ]\ where Hj^ is the Jacobian matrix (e.g., directional cosine matrix), 
and where [Ax, Ay, Az, cAx ]^ are corrections to the baseline coordinates and clock offsets 
of the receivers. Thus, we will model the above equation as: 

(95^ ~ 95^) - A"^ •(/>/ - />/) = iV + A"^ • Hj^ lAx, Ay, Az, cAr f (15) 
The pseudoranges will be used to estimate the vector Hj\Ax, Ay, Az, cAx ]^ as follows: 

ilf-yf)- iPf - Df) = H/iAx, Ay, Az, C'At ]^ (16) 
This will be done through the formation of observation vectors, state vectors, and 
observation matrices, and an estimation process over a plurality of time moments, as 
described below in greater detail. Vector iV will be jointly estimated in this process. 



Generation of the Vector of Observations. 

A vector of observations ^ is generated at each clock time moment tj=j'ATj and 
comprises l-n components, where n is the number of the satellite chaimels. The first n 
components are the residuals of the single differences of the Base and Rover pseudo- 
ranges, which we denote in vector form as A^-: 

^Yi-ri-J>i (17) 

with Yj^Tf-jf, andDy =Df-Df: and 

where ^ and yf are vectors containing the pseudo-ranges of the satellites, 
measured in the Base and Rover receivers, respectively; and 

R B 

Dj and Dj are vectors of the estimated ranges of the satellites fi-om Base and 

R B 

Rover stations at the moment of j-th signal radiation Dj and Dj are computed 

from the known positions of the satellites, the estimated position of the rover, and 
the known or estimated position of the base station by standard geometry. 
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The next n components of the observation vector /Jj are the residuals of the single 
differences of the Base and Rover full phases (A^- ), which we indicate here in vector 
form: 



-1 



Ag>j = (pj-X Dj (18) 
5 where (Pj^ <Pj - (Pj \ 

(Pj and q>j are vectors of the full phases of the given satellite signal, measured at 
the Base and Rover receivers, respectively (the phases are measured in cycles); 

and 

.-1 . 

A IS a diagonal matrix of inverse wavelengths, where each diagonal component 
10 corresponds to a channel and is equal to l/\ where Xis the wavelength of the 

carrier signal in the given channel. 

The fiill phase vectors and^y may be constructed in the form provided by equation 
(5): 

(Tj) -fp.nom\Tj - r/) - (tfj^NCO{Tj\ 
1 5 q^J (Tj) =fp,nomiTj - T/) - (f/j^NC0(Tj), 

or may be constructed in the form which includes the phase offsets of the base receiver: 

(Tj) = \fp.nomiTj - r/) - <ifj,NC0{Tj)\ - <t)'f , 
(fj (Tj) = \fp,nom(Tj - T/) - (lfj,NCO(Tj)] - <l}'/ . 

In either case, we will use the convention practice and have the base receiver correct its 
20 clock so that the base time Tj will be equivalent to the GPS time tj for the purposes of 
implementing our inventions (the time offset has already been accounted for in the above 
equations). For all practical purposes, times Tj and tj will refer to the same processing 
time increment, and can be interchanged in the above equations. 



State Vector Representation. 

The forms represented by Equations (15)-(18) maybe represented in matrix form 
equivalent to: 

IJj = HfAp (19) 
22 



where vector Aj is a state vector related to observation matrix f£f, and where matrix H/^ 
is an observation matrix that specifies the relationship between the components of the 
observations vector and the state vector yi;. State vector vector comprises (4+n) 
components. The first three components are increments (Ax, Ay, Az) to the coordinates 

, JK** , 2:** ) of the baseline vector unknown at the 7-th clock moment, the fourth 
component is the vmknown increment of the reference oscillator phases (c Ar), The 
remaining n components are the unknown floating ambiguities, different in different 
channels (iV', A^^... IsT), Matrix//)^ comprises 2n rows and (4+n) columns, and may be 
divided into the following 4 parts (sub-matrices): 



'nxn 



I 



(20) 



The first part, the left upper comer of this matrix (the first four columns by the first n 
rows), is occupied by the observation matrix relating to the pseudo-range 
measurements, each row corresponding to one channel (fi-om the 1-st to the n-th). For the 
w-th channel, the corresponding row appears like this: 

f f Pjn f hjfi 9 1] , 

where ajn , fijn , hjn - the directional cosines of the range vector to the n-th satellite fi*om 
Rover for the j-th time moment. Methods of computing directional cosines are well 
known to the art and a description thereof herein is not needed for one of ordinary skill in 
the art to make and use the present invention. 

The second part of matrix H/^, the left lower comer (the first four columns by the 
last n rows), is occupied by the matrix product A"^ ^H/^ relating to the full phase 
measurements, each row corresponding to one channel (firom the 1-st to the w-th). For the 
n-th channel, the corresponding row appears like this: 

where Xn is the wavelength of the carrier signal in the n-th channel. 
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The third part of matrix Hj^^ the right upper comer (the last n columns by the first 
n rows) is occupied by zeroes. And the fourth part, the right lower comer (the last n rows 
by the last n columns) is occupied by the elements relating to the floating ambiguities. 
This part represents the identity matrix I^, with dimensions of nxn. For the discussion that 

follows below, it will be convenient to identified sub-blocks of matrix Hf^ as follows 
= [q j I g] , where Qj is a compound matrix formed by and A as follows: 



(having dimensions of 2n x 4), 



and where G is a compound matrix the zero matrix 0„x« and the identity matrix I„ as 
follows: 



G = 



. I« J 



(having dimensions of 2n x « ). 



Equation (18) has 2/2 equations (according to the number of components of the 
observation vector), and may be used to determine the state vector Aj at the y-th clock 
moment (/.e., to determine 4h-« unknown values). Solution of such a system of equations 
at n>A may be performed by the method of least squares. However, our invention relates 
to solving for the ambiguity vector N, which is a component of Aj^ over a plurality of time 
moments. Before we describe that process, however, we want to first describe the 
groundwork of how we can integrate the minimization of the cost fimction F(*) for 
constrained distance between receivers, if used, with the solution of Equation (19), and to 
then describe some covariance matrices that characterize the accuracy of the measured 
data used to generate These covariance matrices are helpful to our invention, but not 
strictly necessary. 

The cost fimction F(*) may be expressed in the following form: 

where: 
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«7 = 



0 = 



Az j 

C-ATjj 



c 

yj 

o 

V 0 J 



is a vector of the first four components of /4y\and 



is a vector of the coordinates of the baseUne vector, with a 



zero as the four vector component. Using a second-order truncation of the Taylor series 
expansion, the cost function F(*) may be approximated as: 



TdF \ T d^F r, 1 



F(r, -ha,)-F(r,) « a/ ^ + - ^^^^ = ""^^^ ^^J^J^J (21) 



where 



ry, and 



1- 



(22) 



(23) 



Because forms (21)-(23) share common variables with equation (18), the minimization of 
10 ^•F(*) can be integrated with the solution of equation (18), as described below in greater 
detail. 
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Generation of the measurements covariance matrix. 

Measurements covariance matrix Rj is preferably formed in the following way on 
the basis of weight coefficients obtained in Base and Rover receivers: 



0 



nxn 



0 



nxn 



(24) 
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where Rj = 



(K^ 
0 

0 

0 
0 



0 

y.2 



0 



0 

(K<P ) 
y,2 



-1 



-1 



0 



0 

0 



(K^ ) 

0 
0 

(K'P ) 



-1 



, and 



~1 



(25) 



The weight coefficients {Kj^ , Kj2 , . . . , Kj^ ) characterize the accuracy of the 

measurements of the residuals A}j of the pseudo-range single differences for the 
5 corresponding satellite channels (1-st, 2-nd, . . ., n-th). The generation of the weight 
coefficients is not a critical part of the present invention, and one may use his particular 
method of weighting. We present here one of our ways, where each of these coefficients 
may be determined according to the weight coefficients measured in each channel by Base 

and Rover receivers for the pseudo-ranges, i.e., by values and A^jf^ . 
1 0 Thus, for example, for the n-th channel: 

where andK^^ are determined taking into accoxmt the measured signal-to-noise ratio 

in the receivers and the satellite elevation angles in the n-th channel (of Base and Rover, 
respectively). Specifically, for each of the receivers (no superscript used), 

15 ^J,m " Zk^m'Sin(^k,rn " ^mird'^y^ when ^k^^ri > ^min , and 

" ^ ^^^^ ^k,^ ^^min, 

where Zfc^^i is the signal strength of the m-th satellite carrier signal as received by the 
receiver (it has been normalized to a maximum value and made dimensionless), where 
(^k,m is the elevation angle of the wj-th satellite as seen by the receiver, where a minimxmi 
20 elevation angle ^min at which the signal becomes visible at the receiver, where Qy^ is the 
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variance of the code measurements (ay»l m). The factor Zitfw*sin(^it,/w ~ ^wm) is 
dimensionless. 

Weight coefficients K^^ , ? . • . 5 AT characterize the accuracy of the 

measurements of the residuals A^y of the phase single differences, and are determined 
similarly. Here the same input data is used: the signal-to-noise ratio and the angle of 



elevation, but another scale for the phase measurements is considered {e.g., instead of 
2 

Qy ). Specifically, for each of the receivers (no superscript used), 

" ^ ^^^^ ^^^^ ^^min, 

where Z^^^, ^k,m and ^min are as they are above, and where (3^ is the variance of the 
code measurements (a(p « 1 mm). 

When the magnitudes of either of weight coefficients K^^^ or K^^^ is less than a 

first selected small threshold value, the value of ^ is generated as a first small 
number which is less than the first threshold value. This is equivalent to setting 



to a large number equal to the inverse of the first small number. Similarly, 



when the magnitudes of either of weight coefficients or is less than a second 

selected small threshold value, the value of is generated as a second small number 

which is less than the second threshold value. This is equivalent to setting {k^ ^ to a 

large number equal to the inverse of the second small number. 

In some instances, satellite signals are blocked from view and should be excluded 
from the ambiguity resolution process. The elements of the covariance matrix Rj 
corresponding to these satellite signals are replaced by a very large number. The very 
large number is selected in advance, and has a value which exceeds by several orders of 

magnitude the nominal values of the covariance matrix components (XJ )~* or (XJ)'^ 
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encountered during operation. Consequently, in further computations, the weights of all 
measurements relating to these channels become so small that they do not influence the 
result. 

As a more simple approach, but currently less preferred, one can use the following 
form of covariance matrix Rj^ or a scaled version thereof: 



1 



0 ^nxn 

'r 



0 



'I 



This form gives an equal weighting of \l<5y of to the rows of pseudorange data provided 

by A^- , and an equal weighting of l/Cq, to the rows of phase data provided hy A <Pj. 
10 When a satellite is blocked from view or not visible, the elements of the covariance 

matrix Rj corresponding to the satellite are preferably replaced by a very large number, as 
described above. 
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Central Aspects of the Present Invention 

Summarizing equation (19), the linearized measurement model at the j-th epoch 
takes the form: 



c- At , 



N 



J J 



f \ 



(26) 



We expect that the vector of floating ambiguities Nj will be constant in time, and 
therefore we will drop the subscript index j. Listead, we will denote Ny the evaluation of 

N obtained through the epoch from 1 to 7 . 
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Generation of the floating ambiguity estimations 



We create the following scalar quantity 7 in equation (27) which integrates the 
minimization of equation (26) with the minimization of the truncated cost functions Fj^ 
for the constrained distance condition (if used) according to equations (21) - (23): 



As mentioned above, q is a scalar weighting factor (penalty parameter) which is equal to 
or greater than 0 (q >0). We seek to minimize / in value in order to obtain 

, A: = 1, . . . , 7 and N , The minimization of scalar J essentially seeks to minimize the 

errors in the individual forms of || Hjt'*Ajt - M^t \f and conditions (1 1), if applicable, for 
data of all of the j epochs being considered, but with the constraint that floating ambiguity 
vector N in each individual state vector A^rbe the same. The individual vectors a^t are 
allowed to have different values. If the fixed distance constraint between rover and base 
stations is not considered to be present, then cost function Fk is omitted fi-om equation 
(27), and all following equations based on equation (27), This may be simply done by 
setting weighting parameter q=0. If the fixed distance constraint is used, values for the 
weighting parameter q can range firom approximately 0.5 to approximately 4, with a 
typical range being between approximately 1 and approximately 3. The best value of q 
often depends upon the amount of noise in the signals and the receivers, and can be foxmd 
by trying several values within the above ranges {Le. , fine timing). The inventors have 
found a value of q=2 to be useful for their test applications. In the case where the 
distance between the receivers is constrained to a fixed value, we emphasize that the use 
of the cost function qF is optional, and that one is not required to use it. Using the 
methods of the present invention without the cost function qF will still provide estimates 
of the floating ambiguities. The inclusion of the cost function qF generally enables more 
accurate estimates. 




(27) 
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The inventors have discovered that a set of N which minimizes scalar Jean be 
found by solving the following block linear system for the set N, which is a vector. 



10 
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Q^R/Ql+gSi 


O 


O 


Q^Rr^G 


O 


QjR2*Q2+gS2 


O 


QjRi^G 


O 


O 








G^R2^Q2 


G^rtIq. 


Xg^r^Ig 

2.-1 



a2 



_i 
N 



T 1 ™ 

A:=l 



To the inventors' knowledge, the form of scalar /provide by equation (27) and form of 
the linear system provided by equation (28) are not found in the prior art. While it is 
preferred to use the weighting matrices R and their inverse matrices, it may be 
appreciated that the present invention can be practiced without them, hi the latter case, 
each weighting matrix R may be replaced by an identity matrix of similar dimension in 
Equation (28) and the following equations; each inverse matrix R'^ is similarly replaced 
by an identity matrix. 

The inventors have further constmcted an inverse matrix for the block matrix on 
the left-hand side of the equation (28), and from this inverse matrix have found a form of 
N which satisfies equation (28) as follows: 



N = 



|](GTRj^lG-G^R^lQ^(QjRi^lQ^+^S^)-lQjR^lG) 



X(GTr^V^ -GTR^lQ^(QfR^lQ^ +^S^)-\QjRi^Vit -^h/t)) 

k=l 



(29) 



where this form comprises annxn inverse matrix multiplying an n x 7 vector. Ih the 
20 discussion that follows, we will identify the w x w inverse matrix as MT^ and the « x 7 
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vector as B, with N = M"^ B. The above form may be more simply expressed if we form 
a matrix matrix Pk for the data of the k-th time moment in the following form: 



and a vector: 



With matrix Pk, equation (29) may now be written as: 

n-1 



N = 



J 



X(G^P;tG) 

k=l 



k=l 



(30) 



We refer to P/^ as a projection-like matrix for q=0 and a quasi-projection matrix for q>0. 



Each component matrix P* G of M= 



is symmetric positive definite 



and may be inverted. Moreover, a matrix comprised of a summation of symmetric 
positive definite matrices is also symmetric positive definite. Thus, matrix M is 
symmetric positive definite and can be inverted. In preferred implementations of the 
present invention, the inverse of M is not directly computed. Instead, a factorization of M 
into a lower triangular matrix L and an upper triangular matrix U is produced as follows: 
M = LU. Several different factorizations are possible, and L and U are not unique for a 
given matrix M. The LU factorization of matrix M enables us to compute the floating 
ambiguities N through a sequence of forward and reverse substitutions. These 
substitutions are well known to the art (cf. any basic text on numerical analysis or matrix 
computations). With symmetric positive definite matrices, one may choose L and U such 
that U=L\ which gives a factorization of M = LL^. This is known as the Cholesky 
factorization, and it generally has low error due to numerical rounding than other 
factorizations. 

With N being generated fi-om this form, the inventors have further found that each 
individual vector can be generated according to the following form: 

ak =(Qk^Rk"^Qk +9S;tr\Qk'^Rk"H^k -GN)-qgk). (31) 

The generation of the individual vectors a^t is optional to the process of generating a set of 
floating ambiguities N. However, when using cost function F(*) and equations (22) and 
(23), one can generate a^ in order to update r^t. The forms of N and a^ which the 
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inventors have discovered enable one to generate a vector N without having to first 
generate the individual vectors a^. 

When using the cost function F(*), Sk and hk (and gk which is derived from hk and 
Sk), an estimate of the baseline vector Tk at time moment "k" has to be generated. As 
5 indicated by equations (22) and (23), both Sk and hk depend upon Lrb, which does not 
normally change, and upon Tk, which can change and often does change. For generating 
Sk and hk, it is usually sufficient to use an estimate of Tk, which we denote at r\J and 
which may be provided by the user or general process that is using the present invention. 
As part of generating the calculated distances Dk and Dk , the user or general process 

10 uses the estimated position of the rover and the known or estimated position of the base 
station. An estimate Tk' can be generated by subtracting the position of the base station 
from the estimated position of the rover. If desired, one can refine the estimate by 
generating ak from equation (31), and then using the combination (fk'+ak) as a more 
refined estimate of Tk in generating Sk and hk- To do this, one may first generate an 

15 estimate N' to N by using equation (30) with Sk=0 and gk=0. Then equation (31) can be 
evaluated using N= Sk=0, and gk=0 to generate refined baseline vectors rk =(rk'+ak), 
and initial values of Sk^ hk'j and gk^ Equation (30) is then again evaluated using the 
initial values Sk', hk% and gk' for Sk, hk, and gk, respectively. The process may be 
repeated again. To speed convergence, 

20 one can consider using prior values of Sk' and gk' in generating ak for k>l as follows: 

ak = (Qk 'Qk + q^k-i (Qk ^ (jik - gn' ) - q^-i ) • 

First Exemplary Set of Method Implementations of the Present Invention 

25 The inventors have discovered a nxmiber of ways to employ the form of N 

provided by Equation (30). For post-processing application, the satellite data may be 
collected for a plurality j of time moments, the matrices and Qjt (and optionally S^) and 
the vectors iik (and optionally h^t and gyt) for each k-th time moment may then be 
computed, and the nxn inverse matrix M"* and the nx 1 vector B may then be generated 

30 and thereafter multiplied together to generate the vector N. The pluraUty of time 
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moments may be spaced apart from one another by equal amoimts of time or by unequal 
amoimts of time. As indicated above, instead of generating the inverse matrix M"^ one 
may generate the LU factorization of M and perfomi the forward and reverse 
substitutions. The factorization-substitution process is faster than generating the inverse 
matrix, and more numerically stable than most matrix inversion processes. If necessary. 
Equations (30) and (31) may be iterated as described above until convergence for 
significantly nonlinear dependency of quantities in equations (22) and (23) on the Rover 
position. 

For real-time applications, IVT^ and B may be initially computed at a first time 
moment from a first set of data (e.g., Ri, Qi, /ii ,R2, Q2, M2) and then recomputed at 
subsequent time moments when additional data becomes available. For instance, we can 
initially compute M and B based on / time moments k=^l to kr=l as follows: 



B/ = 



i:(GTp^G) 

I 

X(G'^Pjfc|i^+^git) 



(31A) 
(31B) 



where the subscript has been used with M/ and B/ to indicate that they are based on the 
/ time moments Af=1 to k=L The floating ambiguity may then be computed from 
N/= [M/~*] • B/, using matrix inversion or LU-factorization and substitution. Data from 
the next time moment /+1 can then be used to generate the updated quantities Mi4-i and 
B/+1 as follows using the previously computed quantities M/ and B/i 



M/+7 = M/ + G^P/^iG 

B/+/ = B/ + G^Yi^x +^g/+l. 



(31C) 
(3 ID) 



With an updated set of floating ambiguities being computed from N/+7 = [Mz+y]"* B/+7, 
using matrix inversion or LU-factorization and substitution. It may be appreciated that 
the above forms of M/+/ and B/+i may be employed recursively (e.g^., iteratively) in time to 
compute updated floating ambiguities N/+7 from previously-computed values of the 
quantities M/ and B/ as satellite data is collected. The recursion may be done with each 
recursive step adding data from one epoch, or with each recursive step adding data from 
multiple epochs, such as provided by the following forms: 

k—m 



M„ = M,+ 



k=l+\ 



(31E) 
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B«, — B/ + 



(3 IF) 



where M/ and B/ are based upon y time moments, and M^m and B^ are based on these time 
moments plus the time moments /+1 through to m, where m > /+1 . While these recursion 
5 methods are preferably applied to real-time applications, it may be appreciated that they 
may be equally used in post processing applications. Furthermore, while one typically 
arranges the time moments such that each time moment h^l occurs after time moment 
it may be appreciated that other ordering of the data may be used, particularly for post 
processing applications. 
10 The steps of the above method are generally illustrated in a flow diagram 40 

shown in FIG. 4. Initially, the data is received (e.g., jf, jf, Df, /)/, q)j^, H/ and 
optionally r/ and Lrb) as indicated at reference number 41. At step 42, the method 
generates, for each time moment j\ a vector A}j of a plurality of range residuals of 
pseudo-range measurements made by the first and second navigation receivers in the form 

15 of: A7j—(j^-T^)'-{Pj^-Df). At step 43, the method generates, for each time 

moment y , a vector of a plurality of phase residuals of full phase measurements made 
by the first and second navigation receivers in the form of: Ag)j = - q)^) — A ^ '{Df^ 
— Dj ), where A is a diagonal matrix comprising the inverse wavelengths of the 
satellites, the set of phase residuals being denoted ^sAq)^, A=l, ..,y. At step 44, the 

20 method generates an LU-factorization of a matrix M or a matrix inverse of matrix M, the 
matrix M being a function of at least A and Hi^, for index k of J?)fc^ covering at least 
two of the time moments j. Exemplary forms of matrix M have been provided above. At 
step 45, the method generates a vector N of estimated floating ambiguities as a function of 
at least the set of range residuals A;^, the set of phase residuals A q>k, and the LU- 

25 factorization of matrix M or the matrix inverse of matrix M. Exemplary forms of vector 
N have been provided above. If the cost function F(*) is included, then steps 44 and 45 
may be reiterated, along with a step of generating Uk for refined estimates of r^, as 
indicated above. 
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It may be appreciated that the above steps may be perfomied other sequences than 
that specifically illustrated in FIG. 4, as long as the data needed to perform a specific step 
is available. For example, parts of steps 41, 42, and 43 may be performed as matrix M is 
being assembled in step 44. 
5 The above method may be carried out on the exemplary apparatus 100 shown in 

FIG. 5. Apparatus 100 comprises a data processor 1 10, an instruction memory 112 and 
data memory 1 14 for data processor 1 10, an optional keyboard/display 1 15 for interfacing 
between data processor 110 and a human user, and a generalized data portal 120 for 

p 5 J? B R Y 

receiving the measured data yf ^ , Dj ^Dj ^q)j ^qjj , Hj, and optionally r/ and Lrb, 

10 each data having been described in detail above. Memories 112 and 114 may be separate, 
or difference sections of the same memory bank. Generalized data portal 120 may take 
any number of conventional forms, such as one or more input/output ports, or one or more 
files stored on disk, tape, non- volatile memory, volatile memory, or other forms of 
computer readable medium. The data can be placed in data portal 120 by any number of 

15 means, such as by a user of the apparatus, or by a more general apparatus that utilizes the 
apparatus of the present invention in carrying out its functions (such as, for example, 
computing precise estimates of the position of the rover or the coordinates of the baseline 
vector). In the former case, the keyboard/display 115 may be used to receive information 
fi-om the user as to when new data has been provided to data portal 120 (this may be 

20 useful in post-processing applications). In the latter case, a more general apparatus may 
comprise the rover station (including receiver channels such as that shown in FIG. 1) and 
radio transceiver for receiving data from the base station. In the case that both the base 
and rover stations are located on the same vessel, the more general apparatus may 
comprise the base and rover stations. 

25 Data processor 110 may be configured to implement the above-described method 

embodiments, such as exemplified by the steps in FIG. 4, by running xmder the direction 
of a computer product program present within instruction memory 112. An exemplary 
computer program product 60 is illustrated in FIG. 6. Computer program product 60 may 
be stored on any computer-readable medium and then loaded into instruction memory 112 

30 as needed. Instruction memory 112 may comprises a non-volatile memory, a 
programmable ROM, and hard-wired ROM, or a volatile memory. 
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Computer program produce 60 comprises five instruction sets. Instruction Set #1 
directs data processor 1 10 to receive the measured data from data portal 120. The 
measured data from portal 120 may be loaded into data memory 1 14 by Instruction Set 
#1 . As another implementation, the data may be loaded into memory 1 14 by subsequent 
5 instruction sets as the data is needed. In the latter case. Instruction Set #1 can take the 
form of a low-level I/O routine that is called by the other instruction sets as needed. As 
such, data portal 120 and data processor 110 under the direction of instruction set #1 
provide means for receiving the measured data for apparatus 100. Instruction Set #2 

directs data processor to generate, for each time moment y, a vector Aj^ of a plurality of 
10 range residuals of pseudo-range measurements made by the first and second navigation 
receivers in the form of: A^- = (j^ - yf) - {Df - Df), As such, data processor 110 
under the direction of instruction set #2 provides means for generating the range residuals 
Ayj for apparatus 100. Instruction Set #3 directs data processor 1 10 to generate, for each 
time moment y, a vector 2< of a plurality of phase residuals of fiill phase measurements 
15 made by the first and second navigation receivers in the form of: A q>j = {q>j - (Pj ) — 

A iPj^ - Df ), where A is a diagonal matrix comprising the inverse wavelengths of 
the satellites. As such, data processor 110 under the direction of instruction set #3 
provides apparatus 100 with means for generating the phase residuals A (py The residuals 
may be stored in data memory 114. 

20 Instruction Set #4 directs data processor 1 10 to generate an LU-factorization of 

. -1 

matnx M or a matnx inverse of matrix M, the matrix M being a fimction of at least A 
and for index k of Hi^ covering at least two of the time moments y . Exemplary 
forms of matrix M have been provided above. Matrix M and its LU-factorization or 
inverse may be stored in data memory 1 14. As such, data processor 110 under the 

25 direction of instruction set #4 provide apparatus 100 with means for generating matrix M 
and its LU-factorization or inverse. Finally, instruction Set #5 directs data processor 110 
to generate a vector N of estimated floating ambiguities as a fimction of at least the set of 
range residuals Ayu, the set of phase residuals A q}k, and the LU-factorization of matrix M 
or the matrix inverse of matrix M. Exemplary forms of vector N have been provided 

30 above. Vector N may be stored in data memory 114. As such, data processor 110 under 
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the direction of instruction set #5 provide apparatus 100 with means for generating an 
vector N, which is an estimate of the floating amhiguities. 

The resulting estimates provided by vector N may be outputted on 
keyboard/display 1 15, or may be provided to the more general process through data portal 
120 or by other transfer means (such as by memory 114 if data processor 1 10 is also used 
to implement the more general process). 



Second Exemplary Set of Method Implementations of the Present Invention 

The inventors have developed additional recursive methods that are generally 
better suited to real-time applications. The second exemplary set of methods generally 
facilitate implementations which require less memory storage space and fewer 
computations. We previously defined a projection-like matrix Pr for the data of the k-th 
time moment as follows: 



Equation (32) was then applied to the form of equation (30) to provide the form: 







-Ir 


J 

1 








1 





k=l 



where M and B are identified in the forms of: 



Mj = 




and B/ = 




>=1 





k=l 



(32) 



(33A) 



(33B) 



The forms of equations (33B) can be expressed in the following recursion form: 

M, = My_i +g'^P,G and By = B^i + G"^ P, jiy + ^g,- (34) 

Then 

Ny = (My_, + G^ P, g)-* (B;_, + G^ P, iij + qgj) (35 A) 

noting that Ny = [M,]"' By for the data from time moments 1 through j, we can write N^i = 
[My-i]"' By_i for the data from time moments 1 through j-\ . The latter can be rearranged 
as My_i Ny_i = Bj-i and used to substitute for the term By-i in equation (35A) to provide: 
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N, = (My.i+G'^P,G) ^ (My_iNy_i +G'Tj\if + q^i) (35B) 

To equation (35B), we now add and subtract the tenn Py G N,_i from the second 
bracketed quantity to obtain: 

Nj = (My-i+G'" Py G)"^ (My_i Ny»i+ G'" Py G N^^i -G"^ Py G Ny^i+G'^ Fj Jly + ^gy) (35C) 

5 Noting that the first two terms of the second bracketed quantity can be factored as 

(My_iH-G^PyG) Nj^i and that the factor (My_i-fG^PyG) is the inverse of the first bracketed 
quantity in equation (35C), it can be seen that the first bracketed quantity multiplied onto 
(My_i+G^PyG) Ny_i is simply Ny_i. Therefore, equation (35C) can be simplified as: 

Ny = Ny_i + (My_i + G^ Py G)"^ (- G^ Py G Ny_i + G'^ Py Jly + qgj) (35D) 

10 The second bracketed quantity of equation (35D) can be fiirther simplified as: 

Ny = Ny_i + (My.i+G'^PyG)-^ (g"^ Py ( JJ,- G Ny_i) + ^g,) 

(35E) 

Using equation (34) this becomes: 
15 N,- = + M/ (g^P,- ( H;- G Ny_,) + q^) (36) 

With this, the following recursive method may be used in a real time application: 

(1) Generate initial values of Mo and No, set epoch counter k to zero (k = Q). 

20 (2) Increment the epoch coimter by one; obtain the data needed to generate the 

matrices and Qjt, and the vector iik- 

(3) Generate R^t and Q^, and fik^ If the constant distance constraint is to be used, 
then also compute Sk, hk^ and gk, and select a value for q greater than zero; otherwise, set 

25 qSk =0 and qgk=0 in steps (4)-(7) below. 

(4) Generate P^t in a form equivalent to: 

P^ =Ri^l - R)^l Qk(Qj R;^^ Qk -^gS^y qJ R^K 

30 (5 A) Generate an LU- factorization of Myt, where M)t is in a form equivalent to 

Mjt = M;t-i G. Exemplary ways of generating LU factorizations are 

described in greater detail below. 
OR - 
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(5B) Generate an inverse matrix M^t of Mjt, where Mjt is in a form equivalent to 
Mjt = M^_i + G. In general, this approach requires more computation 
than the approach of step (5 A), and is currently the less preferred approach. 

(6) Generate ^k^'^k-i + Mjt"* (g^P^ ( [i^- G Na:^i) + ^g)t). If the inverse 
matrix M^"* has been generated according to step (5B), this form may be 
generated by conventional multiplication of M^~^ onto a vector F^, where 

= (g^ Vkiv^k -G Nife_i) + ^g^). If the LU factorization for M^fc has been 
generated according to step (5 A), then the quantity Xk^ M^~^ (g^ P^ 
G Njt_i) + ggk) may be generated from a forward and reverse substitution on 
the form: hkVkXk^ Fjfc* where L.^ is the lower triangular matrix and \Jk is the 
upper triangular matrix of the factorization. The quantity Xk^^ then added to 
N^_i to form . 

(7) Optionally generate ikk= (q/ R*"^ Q*)""^ (Qm^ (lik-GNk) -ggk) 

(8) Reiterate steps (2H7). 

For initial values it is convenient to use Mo = 0 and No = 0. One can also use satellite 
data from a single epoch to generate values for Mo and No using forms (4), (5A), (7 A). 
This is in fact what is done with the above steps (2)-(6) are performed on the data for the 
first epoch with Mo = 0 and N© = 0. In step (3), if the distance constraint is to be used, 
one may use r^' to generate Sk, hk? and gk , or one may generate more refined estimate as 
(rk'+ak), with ak being generated as: 

ak =(Qk'^Rk"^Qkr^Qk'^Rk"^(»*k -GN^„i)),oras 

ak =(Qk^Rk"^Qk ^^Sk_0~\Q]J^k'h\^l. -GN^t-O-^'gjt-i). 

Diuing the first few recursion steps, matrix M^ may be ill-conditioned, and thus 
the generation of the LU factorization or inverse of M^t may incur some rounding errors. 
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This problem can be mitigated by using Ra: = I during the few recursions, or using an 
with more equal weightings between the psuedorange and phase data, and then switching , 
a desired form for 

The steps of the above method are generally illustrated in a flow diagram 70 

p 5 i? B R 'Y 

5 shown in FIG. 7. Initially, the data is received (e.g., (Jj , Dj , Dj ^ Hy and 

optionally rf and Lrb) as indicated at reference nimiber 71. At step 72, the method 
generates, for each time moment y, a vector A J^- of a plurality of range residuals of 
pseudo-range measurements made by the first and second navigation receivers in the form 
of: ^Yj —{yf - }f) - (Pj^ - Df\ At step 73, the method generates, for each time 
10 moment y, a vector 2i of a plurality of phase residuals of full phase measurements made 
by the first and second navigation receivers in the form of: ^4 ^ = ~ q)^) — A (Pf 
- Dj ), where A is a diagonal matrix comprising the inverse wavelengths of the 
satellites, the set of phase residuals being denoted as^d , ^1, 

At step 74, the method generates, for time moment y = 1, an LU-factorization of a 

15 matrix Mi or a matrix inverse of matrix Mi, the matrix Mi being a function of at least A 

1 y 

and Hj \ Any of the forms for M described above may be used. At step 75, the method 
generates, for time moment y = 1, a vector Ni as a function of at least Aj^, A q>u and the 
LU-factorization of matrix Mi or the matrix inverse of matrix Mi. At step 76, the method 
generates, for an additional time moment j ?^1, an LU-factorization of a matrix My or a 

20 matrix inverse of matrix M,, the matrix My being a function of at least A and H^. At 
step 77, the method generates, for an additional time moment j ?sl, a vector Ny as a 
function of at least A^-, A qjj , and the LU-factorization or matrix My or the matrix inverse 
of matrix My. At step 78, the estimated ambiguity vector Ny is reported. If the estimated 
ambiguity vector has not achieved sufficient accuracy, as set by the user, steps 76 and 77 

25 are repeated, with steps 76-78 forming a loop. If the estimated ambiguity vector has 
achieved a sufficient degree of accuracy, or it steps 76-78 have been repeated for a 
maximum number of times set by the user, the process is ended. If the cost fimction F(*) 
is included, then step 76 may include the generation of for refined estimates ofr^, as 
indicated above. 



It may be appreciated that the above steps may be perfomied other sequences than 
that specifically illustrated in FIG. 7, as long as the data needed to perform a specific step 
is available. For example, parts of steps 71-73 may be performed as matrix Mi is being 
assembled in step 74, and as matrix Mj is being assembled in step 76. It may also be 
5 appreciated that the data set provided to the process may span several tens of time 

moments to several thousands of time moments, and that the time moment j=l illustrated 
above may correspond to any of the time moments in the data set, not just the earliest time 
moment or the first received time moment. It may also be appreciated that steps 74-75 
may generated their respective M matrices and N vectors from data measured over two or 

10 more time moments as discussed above with regard to equations (31 A)-(31F), as well as 
just one time moment, and that steps 76-77 may generated their respective M matrices 
and N vectors from data measured over two or more time moments as discussed above 
with regard to equations (3 1 A)-(3 IF), as well as just one time moment. It may be fiirther 
appreciated that approproate ones of the appended claims are intended to cover these 

15 possibilities. 

The above method embodiments may be carried out on the exemplary apparatus 
100 shown in FIG. 5, using a different set of instructions in memory 112. Specifically, 
data processor 110 may be configured to implement the above-described method 
embodiments exemplified by the steps in FIG. 7, by running under the direction of a 

20 computer product program present within instruction memory 112. An exemplary 

computer program product 80 is illustrated in FIG. 8. Computer program product 80 may 
be stored on any computer-readable medimn and then loaded into instruction memory 112 
as needed. Instruction memory 112 may comprises a non-volatile memory, a 
programmable ROM, and hard-wired ROM, or a volatile memory. 

25 Computer program produce 80 comprises eight instruction sets. Instruction Set #1 

directs data processor 1 10 to receive the measured data from data portal 120. The 
measured data from portal 120 may be loaded into data memory 1 14 by Instruction Set 
#1. As another implementation, the data may be loaded into memory 1 14 by subsequent 
instruction sets as the data is needed. Li the latter case. Instruction Set #1 can take the 

30 form of a low-level VO routine that is called by the other instruction sets as needed. As 
such, data portal 120 and data processor 1 10 under the direction of instruction set #1 
provide means for receiving the measured data for apparatus 100. Instruction Set #2 
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directs data processor to generate, for each time moment y , a vector Ajj of a plurality of 
range residuals of pseudo-range measurements made by the first and second navigation 
receivers in the form of: A)^= (j^ - ]^)- {Dj^ - Df). As such, data processor 110 
under the direction of instruction set #2 provides means for generating the range residuals 
5 Ajj for apparatus 100. Instruction Set #3 directs data processor 1 10 to generate, for each 
time moment j\ a vector of a plurality of phase residuals of full phase measurements 
made by the first and second navigation receivers in the form of: A ^j^^tpj - <Pj ) — 

A '{Df - Df ), where A is a diagonal matrix comprising the inverse wavelengths of 
the satellites. As such, data processor 110 under the direction of instruction set #3 
10 provides apparatus 100 with means for generating the phase residuals A q)j. The residuals 
may be stored in data memory 114. 

Instruction Set #4 directs data processor 1 10 to generate, for time moment y = 1, 
an LU-factorization of a matrix Mi or a matrix inverse of matrix Mi, the matrix Mi being 

a fiinction of at least A and Hj^, Exemplary forms of matrix Mi have been provided 
15 above. Matrix Mi and its LU-factorization or inverse may be stored in data memory 114. 
As such, data processor 110 under the direction of instruction set #4 provide apparatus 
100 with means for generating matrix Mi and its LU-factorization or inverse. Instruction 
Set #5 directs the data processor 1 10 to generate, for time moment y = 1, an estimated 

ambiguity vector Ni as a fimction of at least A;^, A and the LU-factorization of matrix 
20 Ml or the matrix inverse of matrix Mi . Exemplary forms of vector Ni have been 

provided above. Vector Ni may be stored in data memory 114. As such, data processor 
110 under the direction of instruction set #5 provide apparatus 100 with means for 
generating an vector Ni, which is an estimate of the floating ambiguities. 

Instruction Set #6 directs data processor 1 10 to generate, for one or more 
25 additional time moments j p^l, an LU-factorization of a matrix My or a matrix inverse of 

matrix M/, the matrix M, being a fiinction of at least A and Hj^. Exemplary forms of 
matrix My have been provided above. Matrix My and its LU-factorization or inverse may 
be stored in data memory 1 14. As such, data processor 110 under the direction of 
instruction set #6 provide apparatus 100 with means for generating matrix Mj and its LU- 
30 factorization or inverse. Because of their similar operations. Instruction Set #6 may share 
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or duplicate portions of Instruction Set #4. Instruction Set #7 directs the data processor 
110 to generate, for one or more additional time moments j t^I, a vector Ny as a function 

of at least A^-, A ^ , and the LU- factorization or matrix My or the matrix inverse of 
matrix My. Exemplary forms of vector Ny have been provided above. Vector Ny may be 
5 stored in data memory 114. As such, data processor 110 under the direction of instruction 
set #7 provide apparatus 1 00 with means for generating an vector Ny, which is an estimate 
of the floating ambiguities. Because of their similar operations. Instruction Set #7 may 
share or duplicate portions of Instruction Set #5. 

Instruction Set #8 directs the data processor 1 10 to report vector Ny as having 

10 estimates of the floating ambiguities, and to repeat Instmction Sets #6 and #7 if vector 
does not have stxfificient (or desired) accuracy, or if it is desired to keep the process going 
even through sufficient accuracy has been reached. The resulting estimates provided by 
vector N may be outputted on keyboard/display 1 1 5, or may be provided to the more 
general process through data portal 120 or by other transfer means (such as by memory 

15 1 14 if data processor 1 10 is also used to implement the more general process). 

Give the detailed description of the present Specification, it is well within the 
ability of one skilled in the GPS art to construct all of the above Instruction Sets without 
xmdue experimentation, and a detailed code listing thereof is not needed for one of 
ordinary skill in the art to make and use the present invention. 

20 

Methods of LU-factorization of Matrix M 
First method. 

We now discuss methods of LU-factorization for matrix M. The factorization 
25 methods may be used in any of the above steps or computer instruction sets where an LU- 
factorization is generated or where an inverse matrix is generated, such as in step (5 A) 
described above, and also used in step (5B) to constmct an inverse matrix for M, although 
such is computationally costly. In one approach of LU-factorization, the matrix Mjt_i 
from the previous iteration is retained for the current iteration, and P^t G fi-om the 
30 current iteration is added to it to form the matrix M^ for the current iteration of step (SB). 
Then, any LU-factorization method may be used to find a lower triangular matrix and 
an upper triangular matrix Vk which satisfies Mjt = La: Ujt . Since Mjt is symmetric 
positive definite, the Cholesky method may be used. This method has good nximerical 
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stability, and generates upper triangular matrix such that it equals the transpose of 
lower triangular matrix L^: Ujt = L^^. With this factorization, Mjt = L/. The standard 
Cholesky method requires a square-root operation for each row of the matrix (w rows 
requires n square-root operations). Such operations may be difficult or time consuming to 
perform on mid-range processor chips. A modification of the Cholesky method 
developed by Wilkinson and Reinsch may be used to avoid these square-root operations. 

hi this method, the factorization of M^t = LDL^ is used, where D is a diagonal matrix. 
We refer readers who are not familiar with this area of the art to the following references 
for further information: 

1. Kendall E. Atkinson, An Introduction to Numerical Analysis, publisher: 
John Wiley & Sons, 1978, pages 450-454; 

2. J. Wilkinson and C. Reinsch, Linear Algebra, Handbook for Automatic 
Computation , Vol. 2, publisher: Springer- Verlag, New York, 1971, 
pages 10-30; 

3. Gene Golub and Charles Van Loan, Matrix Computations, publisher: 
The Johns Hopkins University Press, Baltimore, Maryland, pages 81- 
86. 

The above described main instruction sets that generate LU factorizations of matrix M 
can be constructed to include instructions that direct data processor 1 10 to carry out the 
above forms of Cholesky factorization of matrix M under this first method. The 
combination of these instructions and data processor 110 provides apparatus 100 with 
means for performing the Cholesky factorizations. 
Second method. 

Instead of using the modified Cholesky method, the following method developed 
by the inventors may be used. The inventors currently prefer this method. Given that we 
have generated the previous factorization M^_i = luk-x Ljfc_i^, we generate a factorization 
Tjk T/ for the quantitiy O'^ P^: G as follows: Tk = C'^ P;t G. Later, we will describe 
how this factorization T^^ may be generated. Using Tk Tk = G^ Vk G, the factorization 
Ljfc Ljt^of Mjt may be written as: 

U L/ = U-\ U-i ^ + T/ . (37) 
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It is known in the matrix computation art that the product of two nxn matrices X and Y 
is equal to the sum of the outer products of the columns of these matrices, as specified 
below: 

^ T 

XY''= ZXsYs , (38) 

5 

where {Xi, X2 , . . Xn } are the colimms of matrix X, and where {y i, y2 » • - Yn } 

the columns of matrix Y. The inventors have applied this general knowledge to their 
development of the invention to recognize that 

^ T 

T,T/= Z^k^s^k^s. (39) 
^=1 

T 

10 where {t^j, tk,2 > • • - j tk,n } ^re the columns of matrix Tjt. Each outer product ^t^^^ 

is an nxn matrix of rank one. It is known from the article entitled "Methods for 
Modifying Matrix Factorizations" by P.E. Gill, G.H. Golub, W. Murray, and M.A. 
Saunders (Mathematics of Computation, Vol. 28, No. 126, April 1974, pp. 505-535) that 
when a rank-one matrix of the form z is added to a symmetric positive definite matrix 
15 A, a previously computed Cholesky factorization of matrix A may be modified with a 

relatively few nimiber of computations (much less than the number of computations need 
to generate a factorization of A+z z^). The inventors have further recognized that 

performing n such rank-one modifications on the previous factorization Ljt_i Ljt-i for 

T 

Mjt_i using z z^ = ^k,s^k,s for 5 = 1 to 5 = w would also require less computations and 

20 would have better accuracy that a new factorization for M^. 

Let us illustrate this approach by first defining a set of intermediate factorization 

matrices , L2 , . . L„ , each of which has the same dimensions as L^_i and L^. We 

now go through a sequence of n rank-one modifications which will sequentially generate 

the matrices through L„ , with the last matrix L„ being in a form which is 

25 substantially equal to the desired factorization Ljt. The first rank-one modification starts 

T 

with any of the matrices ^k^s^k^s - Without loss of generality, we will start with s=l in 
order to simplify the presentation. We then write: 
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Li tf = Ljfc_i U_r + ^k,l^k,l (40) 



which can be factored as: 

Li Lf = L*_i( I + Cjcf ) U-i^ , (41) 

where vectors Ci and t^^i are related to one another as follows: l^k-l Ci = tk,i. Vector Ci 
5 is readily obtained from vector tk,i with a forward substitution process with the 
previously computed matrix L^t-i, since L^_i is lower triangular. Then, the above 
reference by Gill, et al teaches how to readily obtain a Cholesky factorization of (I + 

CjCi ), which we denote here as Lj Lj . The reader is referred to that reference, as well 

any other references teaching such factorizations, for the implementation details. From 
10 this it can be seen that the form of equation (35) becomes: 

LlLf=U_i(LilT)L,_i^, (42) 

and thus = Lit-l . The multiplication of two lower triangular matrices is relatively 

easy to perform and computationally inexpensive. 

We how perform the second rank-one modification using rank-one modification 
T 

1 5 the matrix t ^ 2*ifc,2 (s=2) in a similar manner by writing: 

l^2^2= -^tk,2^k,2 (43) 

which can be factored as: 

L2L2=Li(i+C2cJ)l| , (44) 

where vector C2 is generated from tk,2 by forward substitution according to: C2 = tk,2- 

T 

20 A Cholesky factorization of (I + C2C2 ) is then obtained, which we denote here as 

1*2 • Th^s, = Lj L2 . The recursion process continues in this manner until L/i is 
computed. The following recursion sequence can be used: 

(1) Allocate matrices Lq , Lj , , . . ., L„ , and Lj , L2 , . . ., L„ , and vectors 

Ci,C2, Cn. 
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10 



(2) Set an index j =1, and set Lq = l^k-l • 

(3) Repeat the following steps (4) - (7) n times: 

(4) Generate vector Cj from vector tkj by forward substitution according to: 

— T 

(5) Generate matrix Ly as a Cholesky factorization of the quantity (I + ^ j ) . 

(6) Generate matrix as the matrix multipUcation of ^j^i and : 

(7) Increment index j. 

(8) Provide L asLj^. 



Third method. 

Close examination of the second method shows that the following more compact 
recursion sequence may be used: 

(1) Allocate matrices L and L and a vectors C. 

1 5 (2) Set an index j =1 , and set L = Tuk-l • 

(3) Repeat the following steps (4) - (7) n times: 

(4) Generate vector C from tkj by forward substitution according to: L Cy == tkj. 

— T 

(5) Generate matrix L as a Cholesky factorization of the quantity (I + cc ). 

(6) Generate the matrix multiplication product L L , and store the result as L 

20 (^-^ j overwrite the storage location of L with the product L L ). 

(7) Increment index j. 

(8) Provide L as L^. 

The above-described n sequential rank-one modifications require approximately half as 
many operations as direct re-factorization. In addition, low rank modifications to a matrix 
25 factorization are generally more numerically stable than direct re-factorization. 
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The previously-described main instruction sets that generate LU factorizations of 
matrix M can be constructed to include instructions that direct data processor 1 10 to carry 
out the above forms steps of factorizing matrix M imder the above second and third 
methods. Specifically, there would be a first subset of instructions that direct the data 
5 processor to generate an LU-factorization of matrix Mj.j in a form equivalent to Ly_i Ly_i^ 
wherein L/_i is a low-triangular matrix and Ly_i ^ is the transpose of L/_i . In addition, 
there would be a second subset of instructions that direct the data processor to generate a 
factorization of Fj G in a form equivalent to T, T,^ = G^ Py G, where T/ is the 
transpose of Ty (examplary methods for this are described below). There would also be a 

10 third subset of instructions that direct the data processor to generate an LU-factorization 
of matrix M, in a form equivalent to L, L,^ from a plurality n of rank-one modifications of 
matrix L/-1, as described above, each rank-one modification being based on a respective 
column of matrix Ty, where n is the number of rows in matrix My. The combination of 
these subsets of instructions and data processor 110 provides apparatus 100 with means 

1 5 for performing the above tasks. 

Generation of matrix Tu 

The last two factorization methods requires finding T^Tk = G^ G so that the 
column vectors tk,s may be used. One can use the Cholesky factorization method to 

20 generate Tk since G^ P^ G is symmetric positive definite (we call this the first method of 
generating matrix T^). The above second subset of instructions be constructed to include 
fiirther instructions that direct data processor 1 10 to carry out the any form of Cholesky 
factorization of G^P^G.to generate Tjt. The combination of these fiirther instructions and 
data processor 110 provides apparatus 100 with means for generate matrix Ty from a 

25 Cholesky factorization of G^ Py G. However, the inventors have developed the following 
more efficient second method of generation of matrix Tjt. It must be noted that this 
second method of generation of matrix Tk is applicable if the penalty parameter q related 
to the constant distance constraints defined by equation (11), and appearing starting with 
the Equation (27), is set to zero. In other words, the second method of generating matrix 

30 Tk described below is applicable if the constant distant constraints are not used. 
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The second method of generating matrix 

The second method is based on a novel block Householder transformation of the 
matrix Pjt G which the inventors have developed. This method is based on 

constraining the weighting matrices and to the following forms that are based on 
a common weighing matrix W: 

where cTy and are scalar parameters selected by the user, and where ^^qps either 

the LI -band wavelength or the L2-band wavelength of the GPS system. If all the 
wavelengths in A are the same, the above forms reduce to: 



(r^)''=-^W, and (r^)-'=-Lw. 

With W, ay, cr^ , 5 R^ and R^ selected, the following steps are 

employed to generate matrix Tk 

1 . Generate a scalar 6 in a form equivalent to: 6 = — (in the case of 



^2 , .2 2 



dual band measurements b = 



). 



2. Generate a matrix H in a form equivalent to H = W H jt- 

3. Generate a Householder matrix Shh for matrix H . 

4. Generate matrix Tjt in a form equivalent to: 



1- 



"QtPS J 



I4 I 04x(/i-4) 



I 



0(„_4)x4 I I(n-4)x(«-4) 



or in the case of dual band measurements: 
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T = 



All I 0„^„ 
A21 I A22 



where sub-matrixces All, A21, and A22 are as follows: 

-Jl-bAi I4 I 04x(n-4) 



All 



= (w(i))2i 



*HH 



I 

0(n-4)x4 I I(n-4) 



A21 



= (W(2))2 



2S 



HH 



I4 I ^4y<(n-4) 



0(n-4)x4 I 0(n-4)x(n-4) 



, and 



A22 



= (w(2))i 



1-& 



<^4x(ii-4) 



I 



0(n-4)x4 I I(n-4) 



The derivation of this method is explained in Appendix A. 

The above second subset of instructions be constructed to include further 
instructions that direct data processor 110 to carry out the above four general step under 
this second method of generating matrix T^. The combination of these further 
10 instructions and data processor 110 provides apparatus 100 with means for generate 
matrix Tj . 

While the present invention has been particularly described with respect to the 
illustrated embodiments, it will be appreciated that various alterations, modifications and 

15 adaptations may be made based on the present disclosure, and are intended to be within 
the scope of the present invention. While the invention has been described in connection 
with what is presently considered to be the most practical and preferred embodiments, it 
is to be imderstood that the present invention is not limited to the disclosed embodiments 
but, on the contrary, is intended to cover various modifications and equivalent 

20 arrangements included within the scope of the appended cleiims. 
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APPENDIX A - Derivation of the Second Method of Generating Matrix T, 



10 



15 



20 



To simplify the presentation of this derivation, we take the case of where the 
measurements are provided in only one GPS frequency band (e.g., Ll-band). 
Generalization to two-band data and multiband data (GLONASS) is straight-forward 

matter. First, let us relate the range weighting matrix (r'^ )r and the phase weighting 

matrix (r^ )~ to a common weight matrix W and two scalar parameters ay and as 
follows: 

(r^)"'=-Vw. and (R*')r'=-Lw, 

where W is a diagonal positive definite n x n matrix. We now want to look at the block 
factorization of the component (q J R^^ )^ , which is part of 

Pjt = R)^^ - R)^^ Qk(Qj Qk)^^ ' To simplify the presentation, we are 

going to omit the time-moment subscript "A:'' and the pseudorange superscript "V from 
our notation for 



, and simply use the notation 



Q = 



H 



in the following discussion. With ttie above, we can write the 



matrix component (q J R^* Qj^ )~ as follows: 



(qT r-1 = 



[h^ J H^A-l] 



W I o 



nxn 



o 



nxn 



I 



W 



H 



-1 



With the condition that the measurement data is within one signal band, all of the 
wavelengths are the same (for the purposes of determining the floating ambiguity), and 
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the diagonal wavelength matrix A may be replaced by the identity matrix I multiplied by 
the scalar wavelength value X: A = XI. This leads to the simplification: 



(qT r-1 q)-' = 



1 1 



-1 

(h^w uy=cri 



2 . jljl 



(h^w H)r^ . 



The quantity 



^2 . ,2^2 



may be represented a single scalar value b{h = 



5 and the above may be further simplified as: 



(qTr-1q)"^ =C72£,(h^Wh)"\ 



10 



15 



20 



The matrix W is diagonal positive definite so we can write W = W"^ W''^ . Then we 
have 

G^P G = ^ W - ^ W — Hcr^6(H^W h)"^ - H^W^ 



,2 X 



I _ A wi/2j, [aT^Mi^xn^^ jjr ^1/2 



w 



1/2 



Denoting the matrix product W H more simply as H (i.e., W H = H ), the above 



form can be rewritten as: 



W 



1/2 



We now generate reduced forms for H and H . We apply the Householder 
transformation to matrix , which is indicated by Householder matrix Shh: 

H^Shh=V = [v I 04^(„_4)} 

where matrix V is a lower diagonal 4x4 matrix. Matrix S is an n x n matrix, and is an 
orthogonal (so that ShhShh^ = In ) since it is a Householder matrix. Multiplying both 
sides of the above equation by Shh\ and using the fact that ShhShh^ = , we can find 
that: 

= [ V I 04x(n-4) ]Shh = V Shh • 
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Then, by the transpose rule, we can find H = V . We now substitute these reduced 

forms for H and in the prior equation for G^PG, and perform the following 
sequence of substitutions, expansions, and regroupings: 



G^PG=-^W^^^ 



ShhShh-4shhV^(vV^^'vS?ih 



w 



1/2 



i«-4v''(v^)"'v-^ V 



I4 I 04x(„_4) 

0(«-4)x4 I 0(n-4)x(rt-4) 



I4 2 '4 ' 04x(„_4) 



I 



0(n-4)x4 I I(n-4)x(w-4) 



Finally, we obtain 



T = — W^^^ S 



HH 



|l--^l4 I 04x(/i-4) 



0(n-4)x4 I I(«-4)x(/i-4) 



In the case when we use both GPS and GLONASS measurements, take matrices 
(R^)"^ and (R^)"^ in the form 

(Rr)-^=J_W, (R^)r^=^AWA, 

where X is the GPS wave length, then all the above reasoning hold including the formula 

b b 
for T . However the value — is close to 1 (but always less than 1), so 1 — j is a small 

value, the operation of taking a square root does not reduce accuracy. 
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This method of generating the matrix T may be extended to LI and L2 band 
measurements. Let , where and A2 



4 



are the wavelengths of GPS LI and L2 bands, respectively, and W is the weighing matrix, 
as above. Then 



W I o 



o 



nxn 



nxn 



w 



1 



o 



w 



(1) 



nxn 



Matrix Q thus takes the form Q - 



H 
H 

A(2)rH 



, so that 



O 



nxn 



J_W(2) 
«2 



, and 



(Ql'R-'Q)=rHT I h'' I hT(a<'))"' I H'r(A<^>)"'"K 



■^nxn 



O 



nxn 



1 o„^„ 


^nxn 1 


Onxn 




^nxn 1 


^nxn 


1 ^nxn 




^nxn 


1 ^nxn 


^nxn 1 





H 

(a4^h 
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2 1 



1 1 

■+ 



H^WH, 



and 



(h^wh)'V 



= C7jZ)(H^WH)r' 



where b is defined as ^ — 



.2 02 ^2 , _2 o2 2i2 • Note that in this casc the matrfx 



10 



G takes the form G = 


02nx2n 








g'''pg = 


— 1 


^nxn 




0„x„ 1 


1 W(2) 



'<p 



9 w I ^nxn 

I 

O I J-w(2) 

^nxn I o 



(A<'))r'H 



I 



o 



nxn I 2 



(w(i))^ 



o 



nxn 



O 



nxn 
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f 


' 1 i " 








V 


/2 



(h^wh)"^ 



H^W2 j W2 



J 



As in the case of LI measvirements, denote W^H = H , and, implementing Householder 
transformation: 



h''Shh=V = [v I 04,(„_4)J that 



IS 



10 we obtain: 



GPG' = 



•9 



15 



20 

















— S 


V 


^2 



HH 



Onxn I 
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Shh I O 



nxn 



'^nxn 



'HH 



0(n-4)x4 



I4 



0(n-4)x4 



7-I4 I 04x(n-4) I T-I4 |04x(n-4) 



*>HH 



(wO))^ 



I o 



nxn 



"^nxn 



^HH 



(W(2))2 



^HH I 



'nxn 



'nxn 



(w(2))^J 



'HH 



xKx 



Shh 



(wO))^ 



I o 



nxn 



o 



nxn 



(W(2))2 



1 0 Matrix K has the following structure: 



K = 



1-^ 
^ 2 



V 1 y 



I4 I 04x(n-4) 



'(n-4)x4 



i(n-4) 



0(n_4)x4 I 0(n_4)x(n-4) 



4 I 04x(n-4) 



0(n-4)x4 I 0(n-4)x(n-4) 



I4 I 04x(n-4) 



0(n-4)x4 I ^(n-4) 



The Cholesky factorization of matrix K: K = KK^ may be obtained using finite formula; 



15 
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K = 



- — -— I 

0(n-4)x4 I 



'nxn 



Hn-4) 



I4 I 04x(n-4) 



I 



[-J2—J1. 



0(n-4)x4 I 0(n-4)x(n-4) 



I4 I 04x(n-4) 
0(n-4)x4 I ^(n-4) 



Thus, for matrix T, where GPG^ = TT'^ , we may write 



T = 



All I 0„^„ 
A21 I A22 



where sub-matrixces All, A21, and A22 are as follows: 

^\-h/)^ I4 I 04x(n-4) 



All 



= (ww)ii 



>HH 



I 

0(n-4)x4 I I(n-4) 



A21 



= (W(2))2; 



'HH 



I4 I 04x(n-4) 



0(n-4)x4 I 0(n_4)x(n-4) 



, and 



A22 



= (w(2))2! 



*HH 



I 

0(n-4)x4 I I(n-4) 
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